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Abstract 

The difficulties that typically prevent numerical solutions from being 
obtained to finite-energy, two-body, bound-state Bethe-Salpeter equations 
can often be overcome by expanding solutions in terms of basis func- 
tions that obey the boundary conditions. The method discussed here 
for solving the Bethe-Salpeter equation requires only that the equation 
can be Wick rotated and that the two angular variables associated with 
rotations in three-dimensional space can be separated, properties that 
are possessed by many Bethe-Salpeter equations including all two-body, 
bound-state Bethe-Salpeter equations in the ladder approximation. The 
efficacy of the method is demonstrated by calculating finite-energy so- 
lutions to the partially-separated Bethe-Salpeter equation describing the 
Wick-Cutkosky model when the constituents do not have equal masses. 
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1 Introduction 



The Bethe-Salpeter equation is a covariant equation that, in some sense, 
is a relativistic generahzation of the Schrodinger equation ahhough it is devel- 
oped from relativistic quantum field theory rather than from relativistic quan- 
tum mechanics. One particularly noteworthy feature of the equation is that 
interactions are retarded so that there is no action at a distance. While the 
Bethe-Salpeter equation is appropriate for studying properties of relativistic 
bound-state systems, heretofore its use has been limited because, even numeri- 
cally, the two-body, bound-state equation has been exceedingly difficult to solve 
13) ■ For this reason various approximations such as the Blankenbecler-Sugar 
approximation (3) or the instantaneous approximation (Q; 0) are often made 
that reduce the covariant equation in four-dimensional space-time to a more 
tractable, approximately-covariant equation in three dimensions. 

If there are no external fields, the Bethe-Salpeter equation is rotationally 
invariant in three-dimensional space so two angular variables can be separated. 
Furthermore, at least in the ladder approximation, the equation can be Wick 
rotated (analytically continued to Euclidean space) 10), which eliminates the 
singularity in the kernel and makes the equation much easier to solve. When 
a Bethe-Salpeter equation has been partially separated and Wick-rotated, it is 
still an integral or differential equation in two variables. The numerical method 
discussed here offers the possibility of obtaining finite-energy solutions even 
when the equation is not completely separable, which is usually the case, and 
and does not require that the masses of the two bound quanta be equal. 

With several exceptions, solutions to the two-body, bound-state Bethe-Salpeter 
equation have been obtained in the ladder approximation only when the equa- 
tion is completely separable or when the masses of the two bound quanta are 
equal. The Wick-Cutkosky model, which consists of two unequal-mass scalars 
interacting via a massless scalar, is completely separable (i^ i .8; ^9; iioj) . 
and the eigenvalue equation for the coupling constant can be solved numerically 
l|fil: lllHl^ . In the zero-energy limit the Bethe-Salpeter equation is rotationally 
invariant in four-dimensional space-time and is therefore separable. Sometimes 
the completely separated equation has been solved numerically. For example, 
Brennan obtained zero-energy, bound-state solutions for two unequal-mass 
fermions interacting via a massive scalar. In the ladder approximation the 
author (^30) calculated zero-energy solutions for a spin-0 and spin-1/2 con- 
stituent with masses that are not equal and are bound by scalar electrodynamics. 

Even if the equation is not completely separable, finite-energy, two-body, 
bound-state solutions can occasionally be obtained if the masses of the two 
bound quanta are equal. For example, Gammel and Menzel f\3) determined the 
bound-state solutions of two oppositely charged fermions that interact through 
minimal electrodynamics. Schwartz (17) and Nieuwenhuis and Tjon deter- 
mined bound-state solutions for two scalars that interact via a third, massive 
scalar. When the two bound scalars have unequal masses, finite-energy solutions 
were first obtained by Kaufmann ifl^ and later by Seto and Fukui H^), who 
reduce the Bethe-Salpeter equation to an infinite system of integral equations 
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in one variable that are solved numerically. In all cases, before the equations 
are solved, they are Wick-rotated jj^) to eliminate the singularity in the kernel. 

Because the energy appears in more than one place in the Bethe-Salpeter 
equation, a Hamiltonian does not exist, and the equation is an eigenvalue equa- 
tion for the coupling constant instead of for the energy. The equation is solved 
by specifying a value for the energy and then, for the chosen value of energy, 
calculating values of coupling constant that satisfy the equation. Although the 
coupling constants are real in the Lagrangian, there are apparently solutions 
to the Bethe-Salpeter equation with complex values of the coupling constant. 
While the Wick-Cutkosky model (0; 0) has only real values for the coupling 
constant, Kaufmann considered two scalars interacting with a third, mas- 
sive scalar and found complex values. Also, for the same equation Seto and 
Fukui 113) found that "there exists a strong indication that complex eigenval- 
ues appear. . . ." Here attention is restricted to solutions of the Bethe-Salpeter 
equation with real values of the coupling constants, which are the more inter- 
esting physically. 

Numerical solutions to the bound-state, Bethe-Salpeter equation are ob- 
tained in five steps: 

1) The singularity in the kernel is removed by a Wick rotation which 
is always possible in the ladder approximation, and is accomplished by making 
the substitution pq ipQ while rotating the path of integration 90° counter- 
clockwise in the complex pQ-plane. 

2) Two angular variables are separated, which is possible because the Bethe- 
Salpeter equation is rotationally invariant in three-dimensional space provided 
there are no external fields. The resulting equation for the Bethe-Salpeter "wave 
function" ^'(ipo,Ps) is an equation in the two variables po and ps = |p|. In 
the ladder approximation a Wick-rotated, partially-separated Bethe-Salpeter 
equation is of the form 

/oo /"OO 
dgo J dqsV{ipo,Ps,iqo,qs)'^iiqo,qs)- (1-1) 

The above equation actually represents Neq equations, where Neq = 1 if both 
constituent quanta have spin zero and Neq > 1 otherwise. Thus, K{ipo,ps) 
and the kernel V{ipo,Ps,iqo, qs) are both Neq x Neq matrix functions. 

3) Zero-energy solutions are calculated. In the zero-energy limit the Bethe- 
Salpeter equation is invariant under rotations in four-dimensional space-time, 
and is, therefore, completely separable. Zero-energy solutions are expanded in 
terms of basis functions that consist of the product of a set of basis functions 

that depend on the magnitude of the Euclidean four- momentum \p\ = 
{Po+p'^y^^ and hyperspherical harmonics in four-dimensional, Euclidean space- 
time. To obtain solutions, each of the basis functions must (very nearly) 
obey the boundary conditions, which are readily calculated ifl^h. Each basis 
function need not obey the boundary conditions exactly provided that a linear 
combination of the basis functions yields a solution that does. 

4) Finite-energy solutions '^(ipQ,ps) are expanded in terms of a set of basis 
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functions {gj{po,ps)}, 




(1.2) 



Two conditions are imposed on the basis functions: (a) The basis functions 
must (very nearly) obey the boundary conditions, (b) A basis system must be 
chosen that, in the zero-energy hmit, devolves to the basis system that yields 
zero-energy solutions. Knowledge of a basis system that yields zero-energy 
solutions provides guidance in constructing a more general basis system required 
to represent finite-energy solutions. 

5) Finally, the partially separated Bethe-Salpeter equation (1.1) is discretized 
by converting it into a generalized matrix eigenvalue equation for the coupling 
constant. One additional condition is imposed on the generalized matrix eigen- 
value equation: In the zero-energy limit the generalized matrix eigenvalue equa- 
tion that yields finite-energy solutions must devolve to the generalized matrix 
eigenvalue equation that yields zero-energy solutions. Discretization can be 
accomplished, for example, using the Rayleigh-Ritz-Galerkin method l|2ll E3) 
or the method of orthogonal polynomials l|14j) . After expressing the solution 
^{ipo,Ps) in terms of basis functions, both sides of (1.1) are multiplied by 
fiPo^Ps) gi[po,Ps)^ and then integrated over the variables pq and ps. The func- 
tion f{po,Ps) niay be omitted or may be chosen so that that the matrices are 
symmetric or have some other desirable property. The integral equation (1.1) 
has then been converted into a generalized matrix eigenvalue equation 



In the above equation, c is a column vector with the elements Cj that are the 
expansion coefficients for the wave function '^{ipo,ps) in (1.2), and the matrices 
Vh and Vah are Hermitian and anti-Hermitian, respectively. Since the Bethe- 
Salpeter wave function is expressed in terms of Nb basis functions as indicated 
in (1.2), (1.3) is an {Neq x Nb) x {Neq x Nb) matrix equation. 

Because there is no obvious way to force the eigenvalues of (1.3) to be real, in 
general it has been extremely difficult to construct a generalized matrix eigen- 
value equation that yields real values for gig2/4:TT that are solutions to (1.1) A 
sufficient condition for obtaining real eigenvalues of a generalized matrix eigen- 
value equation (1.3) is that Vah = 0. K be Hermitian and either K or Vh be 
positive definite. (See, for example, 1)2311 .') In (1.3), K is often Hermitian. But 
if K is also positive definite, then Vah is usually non-zero, and if Vah is zero, 
then neither K nor Vh is usually positive definite. And even if an eigenvalue 
of (1.3) happens to be real, especially when the basis functions do not obey 
the boundary conditions, the eigenvalue typically is not an eigenvalue of the 
Bethe-Salpeter equation (1.1). 

Solutions to some partially separated Bethe-Salpeter equations have been 
obtained when the masses of the two constituents are equal because, in this 
case, the matrix K in (1.3) is both Hermitian and positive definite, and the 



Kc = 



9192 



{Vh + Vah)c. 



(1.3) 
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matrix Vah vanishes because it is proportional to the difference of the masses 
of the two constituents. For example, for the equal-mass and zero-energy cases, 
the Bethe-Salpeter equation describing the Wick-Cutkosky model (J5; 6) can be 
converted into a matrix equation of the form (1.3) where both K and Vh are 
real, symmetric and positive-definite and Vah = because it is proportional to 
the mass difference of the two bound quanta as well as the energy of the bound 
state I2I. 

When a matrix eigenvalue equation is constructed such that the conditions 
discussed in steps 1) - 5) are satisfied, all eigenvalues usually are not real. But 
real eigenvalues are obtained, and almost all real eigenvalues are solutions of 
the original Bethe-Salpeter equation. 

To demonstrate the techniques for solving a finite-energy, two-body, bound- 
state Bethe-Salpeter equation as well as the effectiveness of the method, finite- 
energy solutions are calculated for the partially separated Wick-Cutkosky model 
when the constituents masses are unequal. Although the equation is separable 
and the solutions were originally calculated from a completely separated equa- 
tion, the method used here only requires that the two angular variables associ- 
ated with rotations in three-dimensional space be separated. The advantage of 
demonstrating the technique with the Wick-Cutkosky model is that the com- 
plications associated with higher spin are avoided. Earlier the author suggested 
an alternative method (0) for solving such equations and then used the method 
to obtain solutions to the the partially separated Wick-Cutkosky model. How- 
ever, the solutions obtained here are more accurate, the numerical method is 
less difficult to implement, is more widely applicable and is more efficient. 



2 The Bethe-Salpeter Equation for the Wick- 
Cutkosky Model 

The Wick-Cutkosky model consists of two scalars with respective masses 

mi and m2 that interact with a third massless scalar. In the ladder approx- 
imation, the Bethe-Salpeter equation that describes a bound state of the two 
massive scalars is 

{(p^ + ^K''){p^ + ^K,yml}{[p-' + {^~ l)if1b. + (e - l)K,] - mjjxKip) 

' -XKiq), (2.1) 



where the notation is that of (25). The parameter < ^ < 1 in the above 
equation is associated with the definition of the center-of-mass variables, and 
K^^ is the four momentum of the bound state. After a Wick rotation (3), in 
the rest frame of the center of mass where K^^ = {E, 0, 0, 0), the Bethe-Salpeter 
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equation takes the form, 

' :XB(«9o,q), (2.2) 



where (p — g) • (p — g) = [pa — qaf + {p^ ~ 9*)(p' ~ <?*) is the Euchdean scalar 
product. 

Dimensionless variables are introduced by defining mi = m(l + A), m2 = 
rn(l — A), dimensionless momentum p' = p/m and dimensionless energy e = 
E /2m. When written in terms of dimensionless parameters, the above equation 
becomes 

{(*p° + 2^ef - p2 - (1 + A)2}{[jp„ + 2(^ - l)ef - p^ - (1 - A)2}xij(*Po, p) 



[p-q)-{p~ q) 



XB(«9o,q)- (2.3) 



where primes have been omitted since all momenta are now dimensionless. 

For compactness of notation, it is convenient to write the coefficient of 
X-e(*POjP) on the left-hand side of (2.3) explicitly in terms of its real and imag- 
inary parts, 

{(*p° + 2^6)2 - p2 - (1 + A)2}{[^po + 2(e - l)ef - p^ - (1 - A)^} 
= DR + iDi. (2.4) 

From (2.4) it immediately follows that 

Db = [pI + p'- ^ee' + (1 + Af][pl + p2 - 4(1 - O'e' + (1 - A)^] 

+ 16^(1 - O^Vo, (2.5a) 

Di = 4epo{-ebo + P' - 4(1 - O'e' + (1 - A)2] 

+ (1 - Oipl + P' - ^ee' + (1 + A)2]}. (2.5b) 

Because Di vanishes both in the zero-energy limit, e = 0, and, if ^ = 1/2, when 
the two constituents have equal masses, A = 0, it is relatively easy to obtain 
solutions in these two limits. 

Since the coupling constant A is real in the Lagrangian, the physically inter- 
esting values of A are real. Actually, for the Wick-Cutkosky model all eigenval- 
ues are real (j^ [hI although for other Bethe-Salpeter equations, solutions 
may exist for complex values of the coupling constant as discussed previously 

Writing Xsiipa, p) in terms of real and imaginary parts, 

X£;(ipo,p) = Xfl.(Po,p) +iX/(Po,p), (2.6) 
and noting that the real and imaginary parts of (2.3) must vanish independently. 
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yields the following two equations: 



n'^m? J_^ {p- q) ■ {p- q)' 



DrXr(j>o,p) - Dixi{po,p) = I — — — -XRiqa^l), (2.7a) 



A f°° d^q 

DiXr{Po,v) + DrXi{Po,v) = I 7 rX/(9o,q)- (2.7b) 

TT^rn^ {p-q)-{p- q) 



Adding (2.7a) and (2.7b), 



DR[xn.{po,p) + XiiPo,p)] + Di[xr{po,p) - X/(Po,p)] 

= ^ L {p-q)\p-q) ^'''^'''''^ + ^'-'^ 

From (2.3) it immediately follows that if XE{ipo,p) is a solution, then 
Xb(~^Po,p) is a solution with the same eigenvalue. Thus, without loss of gen- 
erality it is possible to choose 

Xsiipo, P) = X*E{-iPo, P). (2.9) 

Taking the complex conjugate of (2.9), 

X*EiiPo,p) ^ XE(-ipO:P)- (2.10) 

Therefore, the real and imaginary parts of the solution can be chosen, respec- 
tively, to be even and odd functions of po- 
Defining 

■0(^0, P) = Xij(Po,p) +X/(Po,p), (2.11) 
it immediately follows that 

V'(-Po, P) = Xr{Po, P) - Xi{Po, P). (2.12) 
Consequently (2.8) can be rewritten as 

DrHpo, p) + DM-Po, P) = ^ y_ (p-,).(p-,) ^(^°' ^'-''^ 
which is in a form that is convenient to solve numerically. 

3 Numerical Solutions to the Partially Sepa- 
rated, Wick-Cutkosky Model 

Since (2.13) is manifestly invariant under rotations in three-dimensional space, 
the angular dependence associated with such rotations separates. The wave 
function V'(PO)P) can be written in the form 

i,{po,p) = Fipo,\p\)yi{o,^), (3.1) 
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where Y^{9,(l)) is a spherical harmonic. The integration over the two angular 
variables on the right-hand side of (2.13) can be performed analytically using 
Hecke's theorem |26), and the angular dependence of the solution then sepa- 
rates. Unfortunately, with this approach the remaining kernel is an associated 
Legendre function containing a logarithmic singularity that is difficult to inte- 
grate over numerically l)27|) . Furthermore, the two remaining integrations on 
the right-hand side of (2.13) must be performed numerically. 

An easier method for solving (2.13) is achieved by first rewriting the equation 
in terms of spherical coordinates in four-dimensional, Euclidean space-time |28l 
E): 

p° = \p\ cos 01 pz = \p\ sin 01 cos 02 

= IpI sin 01 sin 02 sin pj, = |p| sin0i sin02 cos0 (3.2) 

The four- momentum is written similarly in terms of primed angles. 

The solution 'i/;(po,p) is then expressed as a series expansion in terms of 
hyperspherical harmonics Plff {cos 0i)y^(02,0) in four-dimensional, Euclidean 
space-time. Defining z = cos 0i, the spherical function P^^^ (z) is given by l(29f ) 

Pi5(z) = (l-zY/^^C^(z), (3.3) 

where Cl{z) is a Gegenbauer polynomial. Now Cl{z) is an even or odd function 
of z if the integer k is even or odd, respectively. From (3.3) it then immediately 

(2) . 

follows that g is an even or odd function of cos 0i if fc — is respectively, an 
even or odd integer. Recalling that xr{P0iP) and xiiPo,p) are, respectively, 
even and odd functions of po, implying that they are also, respectively, even and 
odd functions of cos 0i , zero-energy solutions can be obtained from expansions 
of the form 

Xi?(Po,p) zero— energy — 

J2 9n G,,{\p\)P^^{,{cos 0i)y^(02, (/.), (3.4a) 

ri=l 
Np 

X/(Po,p) zero— energy — 

J2 9n Gn{\p\)P',%+,/cOS 0l)il(02, </-). (3.4b) 
n=l 

In the above expansions, gn is an expansion coefficient, the index i = 0, 2, . . . 
is an even integer and {G„(|p|)} is a set of basis functions, each of which (very 
nearly) obeys the boundary conditions. 

A generalization of the zero-energy basis system (3.4) that is suitable for 
calculating finite-energy solutions is 

Np 

Xr{P0,p) = Y. E 5n,fcG„(H)Pi5(co5 0i)y,l(02,</.), (3.5a) 

n=l k=i,i+2,... 

Np -ff„ax 

X/(P0,P) = E E 9n,kGn{\p\)Pl^}{cose,)Yi{e2A)- (3.5b) 

n=l k=t+l,t+'i,... 
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The values of the index k in (3.5) are chosen so that xr{po, p) and x/(Po, p) are, 
respectively, even and odd functions of cos 9i. In the above expansions, gn,k is an 
expansion coefficient and {G„(|p|)} is a set of basis functions, each of which (very 
nearly) obeys the boundary conditions and will be specified later. Recalling that 
V'(POi p) = Xr(Po^ p) + Xi{po^ p)? the expansion for ip{po, p) immediately follows 
from (3.5), 

V'(Po,p) = E E 3n,fcG„(H)Pi5(cos0i)yj,(^2,</'). (3.6) 

Ti=l k=l,l+l,... 

In the zero-energy limit, the angular dependence of the solution separates 
and only X-rCpojP) or XiiPo^P) is nonzero. Thus, zero-energy solutions can be 
obtained from (3.6) by choosing one value of the parameter k = £,£+!,..., 
with each different value of k yielding different solutions. As a consequence, in 
the zero-energy limit the basis system (3.6) devolves to a suitable basis system 
for obtaining zero-energy solutions. 

There are three advantages to seeking solutions of the form (3.6) instead 
of (3.1): 1) After using Hecke's theorem |26l) to perform the three angular 
integrations analytically, the remaining kernel does not contain a logarithmic 
singularity. 2) In (2.13) only one integration must be performed numerically 
instead of two. 3) The basis functions have the correct angular dependence 
for zero-energy solutions so that fewer angular terms are required to obtain 
accurate, finite-energy solutions when the states are tightly bound. 

Substituting (3.6) into (2.13), 



E E 9n,kGn{\p\)[DRP\^}{coser) + DiP^^}{-cos9^)]Yi{e2, 

n=l k=tj+l.,... 

, Np K^^^ oo 

^E E / dM|gpG„(M) 

n=l k=e,i+i,... •'^ 



{p^ + — 2pqcos Q) 



^A'hcos' e,)Yi{e'^A% (3.7) 



where Q is the angle between the four-vectors p and q. Using Hecke's theo- 
rem |2^ to perform the angular integration (All necessary formulas are in the 
appendix of Ref. (,29.1 .). 

J2 E 5n,fcG„(b|)[i?flP£)(cos 0i) + DjPI:'}{-cos 9,)]Yiie2, <j>) 

n=l k=i,£+l,... 

. Np Krr... ^OO 

= ^E E 9n,k d\q\\q\'G^{\q\)Ai'\\p\,\q\)Pl:'}{cose,)Yi{e,,^ 

7T Tfl In 

n=i k=e.e+i.... •'^ 

(3.8) 
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The function A^^^(|p|, is (|2^|2gD, 



(3.9) 



The dependence on the angular variables 02 and (j) separates as it must. 

To determine the boundary conditions, the parameters go and goo must be 
calculated that satisfy 



G„(k|)^M*% 
G„(kl) \q\-'-. 
Once the asymptotic behavior of integrals of the form 



Hp) 



d\q\\qrG{\q\)Ai'\\pl\q\) 



(3.10a) 
(3.10b) 

(3.11) 



which appears in (3.8), are determined, the boundary conditions are readily 
calculated. Specifically, the parameters iq and «oo must first be calculated that, 
respectively, satisfy 



i{\p\)-^\pr, 

|p|->0 



im — \p\- 



(3.12a) 
(3.12b) 

There are two possible values for the parameter iq in f3.12a') l(l5|) : 

Solution lA : Iq — k, — n + k + 1 < go (3.13a) 

n - fc - 1 < goo 

Solution IB : iq = go + n — 1, — n — k— l<go< —n + fc + 1 (3.13b) 

n - fc - 1 < goo 

Similarly, there are two possible values for the parameter ioo in (3.12b) ifisl) : 

fc + 2, -n-k-Kgo (3.14a) 



Solution IIA 
Solution IIB 



n + l, 



— n — k — 1 < go 
n + fc + 1 < goo 

-n-k-Kga (3.14b) 
n — fc — 1 < goo < n + k + 1 



Using the fact that as \p\ — > 0, Dji — > constant, Dj \p\ and substituting 
(3.10a) into (3.8), for Solution lA it follows that 



(3.15) 
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Thus, go = k. Because the smallest value of k=£, 



G„(b|) ^ \pf. (3.16) 

IpHo 

As can be readily checked, there are no other solutions for Qq. Similarly, at large 
\p\ , the only solution is ^oo = + 6 so 

^"(l^'D.^pTe- (3-17) 

The knot structure is as follows: There are Np cubic splines in the expansion 
(3.6) and Np + 4 momentum knots Tp{i). To determine the momentum knots, 
Np Chebyshev points Xp{i) are calculated on the interval —1 < Xp{i) < 1, 

Xp(i) = -cos^2^^^, i = 1,2,..., Np. (3.18) 

The momentum knots Tp{i + 4) are then given by 



Tp{i + 4) = C'^^±^^+C", I = 1,2,..., Np. (3.19) 

The constant C is chosen by trial and error to approximately minimize the 

lowest zcro-cncrgy eigenvalue, and the constant C" is chosen so that the first 
knot on the positive |p|-axis is not too close to \p\ =0. The values C = 1.0 and 
C" = 0.01 were satisfactory. A knot is placed at the origin, Tp(4) = 0, and the 
three knots Tp(l), Tp{2) and Tp{3) are placed on the "negative" \p\ axis to allow 
maximum freedom in constructing the solution from splines near \p\ = 0. The 
three knots on the "negative" \p\ axis are mirror images (about the origin) of 
the first three knots in (3.19). With this choice of knots, the first three splines 
are finite at the origin, creating sufficient freedom to construct solutions from 
splines near |p| = 

Angular knots are chosen on the z axis, where z = cos^i, so that numerical 
integrations can be carried out over cos^i. Defining Ng = Kmax — £+^, which is 
the number of hyperspherical harmonics in the expansion (3.6) of the solution, 
arbitrarily, but in analogy with splines, the number of angular knots is chosen 
to be Ng + 4. The angular knots Tz{l) = -1, Tz{Ne+4) = 1 and the remaining 
knots are the Chebyshev points 

Teii +l)=-cos-^^^-^y i = 1,2,..., No + 2. (3.20) 

So that the basis functions G„(|p|) asymptotically vanish as indicated in (3.16) 
and (3.17), G„(|p|) is chosen as follows: 



GnM = J^Bndpl) ^ GeMBM, (3.21) 
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where ^^Hpl) is a convergence function, "a" is a constant and Bn{\p\) is a cubic 
spline (|30). At small \p\, G„(|p|) ~ |p|^i3„(|p|). Since the splines are also 
functions of \p\, at small \p\ the individual basis functions Gndp]) do not exactly 
obey the boundary condition (3.16). However, as \p\ 0, for each solution the 
sum of the first three splines in the expansion (3.6) approaches a constant so 
that each solution satisfies the boundary condition exactly. At large \p\, since all 
splines vanish, the convergence function is chosen to vanish as l/jpl^"*"^, which is 
one power of \p\ slower than the rate in (3.17). But at large \p\, because the last 
spline does not decrease exactly as l/|p|, basis functions very nearly, but do not 
exactly, satisfy the boundary condition. Solutions can be obtained that obey the 
boundary condition exactly at large \p\ by extending momentum knots beyond 
\p\ = oo (15) just as they were satisfied exactly by extending momentum knots 
below \p\ — 0. Because solutions decrease so rapidly at large momenta, the 
value of solutions at very large \p\ has minimal impact on numerical solutions. 
As a consequence, for a given number of splines, more accurate solutions are 
obtained without using a knot structure that extends beyond |p| = oo and has 
fewer splines at small \p\ where the solution has most of its support. 

To solve (3.8), the dependence on 62 and (j) is first separated. Then the re- 
sulting equation is discretized using a hybrid method: The angular dependence 
is discretized using the method of orthogonal polynomials (jl41 . which requires 
that the coefficient vanish independently for each of the first Ng spherical func- 
tions Pf^j^__^ ^{z), Ig — 1,. . . ,Ng in the equation. The product of functions 
that appear in the equation and spherical functions that appear in the expan- 
sion for solutions can be reexpressed as spherical functions, some of which have 
a larger first index. As a consequence, although there are Ng different spheri- 
cal functions in the expansion for the solution, there are more than Ng different 
spherical functions in the equation. Consequently, if a solution is to be obtained, 
the series must converge . Using the orthogonality relationship for the spherical 
functions Pif, 

/_;.,(i-=')*p,'->(=)p,;->(=) . (3.22) 

it follows that multiplying the equation by \/l ~ z'^ ^i+ie i ^^"^ integrating 
over z achieves the desired discretization. The momentum dependence is dis- 
cretized using a modified Rayleigh-Ritz-Galerkin method (21; ,22,). Thus (3.8) is 
converted into a generalized matrix eigenvalue equation of the form Ag = -^Bg, 
where the elements of the column vector g are the expansion coefficients gn_k in 
(3.6), by multiplying (3.8) by 

\pfVl^Ql{\p\)BIM)Pwe-^M^ (3.23) 
and integrating over z and \p\. 
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X 

and 



The expressions for the matrices Aij and Bij are, respectively, 

/I poo 
^dzVT^ d\p\\pfg,{\p\)BiM)Pwe-iM 

[DrP^Z.-.A^] + DiP^l\^^,/-z)]BjMSd\p\), (3.24a) 



d\q\\q{- 



/I poo pa 

dzVl-z' / dbl / 
-1 JQ Jo 



(3.24b) 



In (3.24b) 



fel if kl<bl 

Ri\p\,\q\)^{ (3.25) 

g if H<kl- 

As compared with (3.6), indices have been changed in (3.23) and (3.24) so 

(2) 

that terms are automatically excluded when j < j in Pij - Here /p = 1, . . . , Np] 
le = 1, ■ ■ ■ , Ng and the index i is given by i = Np{Ig — l)+Ip with a corresponding 
expression for j. 

With the aid of the orthogonality relationship (3.22) for spherical functions, 
the integral over the variable z in (3.24b) can be performed analytically yielding 







R(\v\ \q\Y+^" 

X g,{\p\)B,M T BjM}Si{\q\) Si, J,. (3.26) 

As can be seen from (3.26), if JV — 3 the matrix B is both symmetric 
and positive definite. Also, the matrix A is symmetric when the quantity Dj 
vanishes, which, from (2.5b), occurs either when the energy is zero or when the 
masses of the constituents are equal. But when A and B are both symmetric and 
at least one is positive definite, all eigenvalues are real ifisj) . so all eigenvalues 
are real for the two cases just mentioned. However, if the energy is finite and 
the masses are unequal, all eigenvalues of the discretized equation are not real, 
but real eigenvalues are obtained. For the solutions of the discretized equation 
corresponding to the lowest six real eigenvalues, which were the only solutions 
checked, when a sufficient number of basis functions were used, the solutions 
to the discretized equation also satisfied the partially separated Bethe-Salpeter 
equation. 

The disadvantage of choosing Af — 3is that the generalized matrix eigenvalue 
equation actually represents the partially separated equation multiplied by |pp. 
The factor of |pp reduces the sensitivity of the matrix equation to the form of 
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the solutions at small \p\ with the result that numerical solutions do not satisfy 
the partially separated equation as well at small Choosing M = 1 allows 
accurate solutions to be calculated for small \p\. But then, even when the energy 
is zero or when the masses of the constituents are equal, all eigenvalues are not 
real. Nevertheless, the same set of solutions is obtained when A/" = 1 as when 
Af — 3. Because the solutions are more accurate at small \p\ when Af = 1, all 
solutions in Tables 1 and 2 are calculated with this value. 

The solutions obtained from the generalized matrix eigenvalue equation 
are checked in two ways: (1) As the number of basis functions 

Gn{\p\) and ^ {cosOi) are increased, the value of each eigenvalue must con- 
verge. (2) For each solution the left- and right-hand sides of (3.7) are compared 
at the center of each rectangle in the physical region of the grid formed by the 
angular knots and the momentum knots. By examining where the left- and 
right-hand sides of the equation agree least well, deficiencies are revealed and 
possible remedies can be efficiently tested. In addition, a reliability coefficient 
Tihs-rhs (31), which is a statistical measure of how closely the left- and right- 
hand sides agree at the {Np) x [Ng + 3) points, is calculated. If the left- and 
right-hand sides agree exactly at every point, then rihs-rhs = 1- 

Table 1 lists values of the coupling constant \/m? that are calculated in the 
zero-energy limit (e = 0) when mi = 4m2. Since the angular dependence sep- 
arates in the zero-energy limit, only one angular basis function Pj^ ^{cos6i) is 
used {No = 1). That single angular basis function is indicated by the value of the 
index fc = ^ in the sum (3.6). As the number Np of momentum basis functions 
G„(|p|) is increased, the calculated values of the coupling constants converge to 
correct values, and the reliability coefficients rihs~rhs approach unity. (The "ex- 
act" eigenvalues are correct to at least four significant figures and are calculated 
numerically from the completely the separated equation (6).) 
Table 1. Calculated values for the coupling constant X/m? in the zero-energy 
limit when toi = 4m2. 





Np = 5 


Np = 10 


Np ^ 20 


-^/'^Gxact 


i 




'^Ihs — rhs 


•^/™calc 


'^Ihs — rhs 


^/^"calc 


'^Ihs — rhs 


1.838 





1.905 


0.9994 


1.841 


0.999990 


1.838 


0.99999968 


5.000 





5.775 


0.9990 


5.035 


0.999993 


5.000 


0.99999972 


5.654 


1 


5.753 


0.9989 


5.647 


0.999996 


5.652 


0.99999995 


9.817 





11.53 


0.9993 


9.996 


0.999993 


9.822 


0.99999957 


10.43 


1 


11.71 


0.9968 


10.42 


0.999991 


10.42 


0.99999989 


11.46 


2 


11.51 


0.9988 


11.43 


0.999986 


11.45 


0.99999982 



Table 2 lists values of the coupling constant A/m^ that are calculated for four 
values of the square of the normalized energy = [E/ {7711+1712)]'^ = 0.1, 0.5, 0.9 
and 0.99 when mi = 4m2. For each energy, the number Np of momentum basis 
functions and the number Ng of angular basis functions used in the calculation 
are listed. As can be seen from Table 2, as the normalized energy e increases from 
zero to unity (and the binding energy decreases to zero), even using additional 
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basis functions it becomes increasingly difficult to obtain accurate eigenvalues. 
Nevertheless, when = 0.99, the first eigenvalue is readily calculated with a 
relative error of a few tenths of a percent, and the first six eigenvalues are all 
determined with relative errors less than five percent. 

Table 2. Calculated values for the coupling constant Xjrv? when the energy is 
finite and m\ = 4m2. 



0.1 Np = 20 Ne = 10 


''^/™cxact 


^/™calc 


Tlhs—rhs 


1.686 


1.686 


0.99999973 


4.690 


4.691 


0.99999975 


5.156 


5.154 


0.99999987 


9.252 


9.264 


0.99999954 


9.688 


9.669 


0.99999960 


10.42 


10.41 


0.9999973 



= 0.5 Np = 20 iVe = 10 


-^/'^cxact 


^/"^calc 


Tlhs—rhs 


1.052 


1.052 


0.9999985 


3.112 


3.111 


0.9999966 


3.344 


3.341 


0.9999963 


6.174 


6.185 


0.9999880 


6.532 


6.493 


0.9999865 


6.748 


6.757 


0.9999858 



= 0.9 Np = 25 iVg = 20 


-^/™cxact 


^/"^calc 


f"lhs~rhs 


0.3167 


0.3165 


0.99985 


0.8500 


0.8487 


0.99973 


1.550 


1.547 


0.99976 


1.590 


1.586 


0.99942 


2.534 


2.522 


0.99899 


2.604 


2.595 


0.99918 



= 0.99 Np = 30 TVg = 30 


'^/™corrcct 


^/"^calc 


'^Ihs — rhs 


0.0702 


0.0700 


0.968 


0.166 


0.164 


0.968 


0.286 


0.273 


0.954 


0.427 


0.415 


0.928 


0.590 


0.613 


0.871 


0.734 


0.718 


0.988 



4 Conclusions 

A systematic method is discussed for solving finite-energy, two-body, bound- 
state Bethe-Salpeter equations that does not require that the equation be com- 
pletely separated or that the constituents have equal masses. To apply the 
method, an equation must first be Wick-rotated |5j) and then the two angular 
variables associated with rotations in three-dimensional space must be sepa- 
rated, which is possible for many two-body, bound-state Bethe-Salpeter equa- 
tions, including all such equations in the ladder approximation. Zero-energy 
solutions are calculated first: The zero-energy equation is completely separated 
by expressing the solution as a product of a hyperspherical harmonic and a 
function that depends only on the magnitude |p| = (p§ -f p^)^/^ of the 

Euclidean four-momentum. The zero-energy solutions are then calculated by 
first expanding the function -Fdpl) in terms of basis functions that (very nearly) 
obey the boundary conditions and discretizing the equation by converting it 
into a generalized matrix eigenvalue equation that is solved numerically. It is 
important to calculate zero-energy solutions first because the basis functions 
that yield zero-energy solutions provide a guide for determining the basis func- 
tions that yield finite-energy solutions. Finite-energy solutions are calculated 
by expanding solutions in terms of basis functions, each of which is a product 
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of a "convergence function" that typically obeys the boundary conditions, a 
hyperspherical harmonic in four-dimensional, Euclidean space and a spline that 
depends on the magnitude of the four-dimensional, Euclidean momentum. The 
basis functions that yield finite-energy solutions must devolve to the basis func- 
tions that yield zero-energy solutions in the zero-energy limit. The partially sep- 
arated equation is then discretized and solved numerically by converting it into 
a generalized matrix eigenvalue equation. The generalized matrix eigenvalue 
equation that yields zero-energy solutions provides guidance in formulating a 
generalized matrix eigenvalue equation that yields finite-energy solutions, and 
the latter must devolve to the former in the zero-energy limit. Even though the 
coupling constants, which are calculated as eigenvalues of the generalized matrix 
eigenvalue equation, usually cannot all be forced to be real, real eigenvalues and 
corresponding solutions are obtained that satisfy the Bethe-Salpeter equation. 

To demonstrate the techniques and utility of the method, when the con- 
stituents have unequal masses, finite- and zero-energy solutions are calculated to 
the partially separated Bethe-Salpeter equation describing the Wick-Cutkosky 
model (ISiIQ). For this particular equation it is convenient, but not essential, to 
discretize the angular dependence using the method of orthogonal polynomials 
()l4|) and the momentum dependence using a modified Rayleigh-Ritz-Galerkin 
method |2lt The advantage of demonstrating the techniques by solving 

the Wick-Cutkosky model is that complications associated with higher spin are 
avoided. 

Using the numerical techniques presented in the paper, the author has begun 
obtaining finite-energy solutions to the scalar electrodynamics model ll^l : Is^l 
and to the scalar-scalar model itl^ list Efil ^ when the bound constituents have 
either equal or unequal masses. Thus, it is highly likely that the numerical 
method discussed here provides a means for obtaining general, finite-energy 
solutions to many two-body, bound-state Bethe-Salpeter equations. 
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